Introduction

In this project, I will be investigating and forecasting electricity consumption in the UK. I will be comparing two different, both very popular forecasting methodologies: Holt-Winters exponential smoothing and ARIMA models. I will be comparing and evaluating both methodologies against the naive model which simply forecasts each observation to be equal to the previous one (Y_t = Y_t-1). I will be using four different performance metrics:

Performance metrics

1.) Mean Squared error (MSE): MSE simply takes the difference between the forecasted value Yhat and the actual value Y, squares it for each observation and then takes the average over all observations. The reason for this is so that negative and positive values don’t cancel each other out. Additionally, this also penalises predictions that are very far off more than ones that are only slighty wrong. 2.) Mean Absolute error (MAE): MAE also computes the difference between forecasted value and actual value, but instead of squaring it simply takes the absolute value of the difference. It then computes the average among all the differences. This method is preferable if we want all errors to be weighted equally, and don’t want errors that are further off to be penalised higher.

3.) Root Mean Squared Error (RMSE): RMSE is simply the square root of the MSE. The reason for this is so the error is on the same scale as the original data, giving better insight into how the error term compares to the data.

4.) Mean Absolute Percentage Error (MAPE): MAPE is the average ratio between the absolute error and the absolute value of the actual value (Y).

Reading in the data

Let’s take a look at the data:

data(USgas)
ts_info(USgas)
##  The USgas series is a ts object with 1 variable and 238 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2019 10
ts_plot(USgas, Xtitle = "Year", 
        Ytitle = "Natural gas consumption in Million cubic feet",
        title = "Monthly US Natural gas consumption")

Looks like there is quite a strong seasonal component in the data, which makes sense as a lot of heating is powered using natural gas so there would naturally be more demand for gas during the winter months compared to the summer months. We can also see a slightly increasing trend from about 2010 onwards. Let’s take a look at a decomposition of the time series and see if this confirms these first visual impressions:

decompose(USgas) %>% plot()

Decomposition reveals a trend as well, beginning around 2005 but becoming more pronounced from about 2010 onwards. I’ll be splitting the data into training and test set, the test set will be the last 12 months of the data. I’ll start off by fitting a naive model, which simply uses the average of all preceding observations. This will be a good benchmark to compare the performance of more complex models to.

USgas_par <- ts_split(USgas, 12)
train <- USgas_par$train
test <- USgas_par$test
naive_mod <- naive(train,h=12)

Fitting a naive model

fc_naive <- forecast(naive_mod,h=12)
accuracy(fc_naive, test)
##                      ME     RMSE      MAE        MPE     MAPE     MASE
## Training set  -1.028444 285.6607 228.5084 -0.9218463 10.97123 1.984522
## Test set     301.891667 499.6914 379.1417  9.6798015 13.28187 3.292723
##                   ACF1 Theil's U
## Training set 0.3761105        NA
## Test set     0.7002486  1.499679
test_forecast(actual = USgas,
forecast.obj = fc_naive,
test = test)

As we can see, the naive model simply takes the average of all previous observations without any considerations of seasonal or trend data. Let’s fit a Holt Winters exponential smoothing model instead and see how well this captures the data and whether we get an improvement over the naive model:

Fitting a Holt Winters exponential smoothing model

hw_mod <- HoltWinters(train)
hw_mod
## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = train)
## 
## Smoothing parameters:
##  alpha: 0.371213
##  beta : 0
##  gamma: 0.4422456
## 
## Coefficients:
##             [,1]
## a   2491.9930104
## b     -0.1287005
## s1    32.0972651
## s2   597.1088003
## s3   834.9628439
## s4   353.8593860
## s5   318.1927058
## s6  -173.0721496
## s7  -346.6229990
## s8  -329.7169608
## s9  -112.1664217
## s10 -140.3186476
## s11 -324.5343787
## s12 -243.9334551

We can see that the model takes into account the average level of the preceding observations (alpha = 0.37) as well as a strong seasonal component (gamma = 0.44). It does not, however, take into account a trend (beta = 0).

fc_hw <- forecast(hw_mod, h = 12)
accuracy(fc_hw, test)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  7.828361 115.2062 87.67420 0.2991952 4.248131 0.7614222
## Test set     51.013877 115.1555 98.06531 1.7994297 3.766099 0.8516656
##                     ACF1 Theil's U
## Training set  0.21911103        NA
## Test set     -0.01991923 0.3652142

This is a significant improvement over the naive model, with RMSE reducing by more than 75% from 500 to only 115. We can see a similar reduction when looking at MAPE and MAE.

test_forecast(actual = USgas,
forecast.obj = fc_hw,
test = test)

Looking visually at our predictions, this confirms what the error metrics tell us. The model beautifully captures the seasonal patterns.

Fitting an ARIMA model

Let’s see how a different class of models, the ARIMA model, performs for this dataset. Arima models rely on a combination of two modelling techniques, the AR and MA processes. The AR process basically explains a model’s future values as a linear combination of it’s previous values, or lags. It requires the time series to be stationary, meaning it can have an increasing or decreasing trend or variance. Lots of time series are not stationary however, be it due to trends, varying seasonal patterns or a certain number of random events. Fortunately, there is a way to deal with these time series as well using a technique called differencing. This mean subtracting the value from one cycle prior (so for example in a daily series, Yt = Yt - Yt-365). After doing this, we simply analyse the difference between these two time points. Forecasts made on these differences easily can be converted back to an actual time series.

Due to the strong seasonality in the data, a SARIMA model, which has three additional parameters to capture seasonality, will likely be the most appropriate. To get a better understanding of the serial correlation of the data (how much each observation is correlated with previous observations, or lags), I will plot the ACF and PACF plots.

par(mfrow=c(1,2))
acf(train, lag.max = 60)
pacf(train, lag.max = 60)

The ACF plot indicates that there is a strong positive correlation with previous observations of the same season, and a weaker but still significant negative correlation with observations of the opposite season. All of this does not come as a surprise since we know that natural gas consumption is highly seasonal. We can also see that the correlation is decaying over time, indicating that the series is in fact not stationary and we will need to do some differencing to continue with our ARIMA model.

USgas_d12 <- diff(train, 12) #specifying the period of differencing, in this case 12 months
ts_plot(USgas_d12,
title = "US Monthly Natural Gas consumption - First Seasonal
Difference",
Ytitle = "Billion Cubic Feet (First Difference)",
Xtitle = "Year")

This yields the following. We have successfully removed the trend in the series, however there still exists non-constant variance which we will need to account for by differencing once more:

USgas_d12_1 <- diff(diff(USgas_d12, 1))
ts_plot(USgas_d12_1,
title = "US Monthly Natural Gas consumption - First Seasonal and
Non-Seasonal Differencing",
Ytitle = "Billion Cubic Feet (Difference)",
Xtitle = "Year")

This looks better, let’s look again at the ACF and PACF plots after the transformation:

par(mfrow=c(1,2))
acf(USgas_d12_1, lag.max = 60)
pacf(USgas_d12_1, lag.max = 60)

As we can see in the plot, after the transformations the correlation with the lags seem to be tailing off very quickly. Let’s try fit an arima model using auto.arima. This automatically determines the best parameters for the AR and MA processes, as well as doing the required differencing for us.

USgas_arima_mod <- auto.arima(train)
USgas_arima_mod
## Series: train 
## ARIMA(2,1,1)(2,1,1)[12] 
## 
## Coefficients:
##          ar1      ar2      ma1    sar1     sar2     sma1
##       0.4301  -0.0372  -0.9098  0.0118  -0.2673  -0.7431
## s.e.  0.0795   0.0755   0.0452  0.0886   0.0826   0.0748
## 
## sigma^2 estimated as 10446:  log likelihood=-1292.83
## AIC=2599.67   AICc=2600.22   BIC=2623.2
USgas_test_fc <- forecast(USgas_arima_mod, h = 12)
accuracy(USgas_test_fc, test)
##                     ME      RMSE      MAE       MPE     MAPE      MASE
## Training set  5.846395  97.81607 73.42635 0.1171733 3.522336 0.6376842
## Test set     37.858513 103.23333 81.47035 1.3112392 3.261787 0.7075437
##                      ACF1 Theil's U
## Training set -0.004151171        NA
## Test set     -0.046726335 0.3404228

Not bad, we reduced RMSE from 115 using HoltWinters to just 103. Let’s take a look at the plot of the forecast:

test_forecast(USgas,
forecast.obj = USgas_test_fc,
test = test)

This model also captures the seasonal component very well, while also including the trend. We can see that the peak for the forecasted time points is just a little bit higher than in the previous years, which is what we would expect looking at previous year’s changes. This is something the Holt Winters model did not pick up on, which is reflected in the slightly lower RMSE. To check how well the model fits our data, let’s perform a residuals check:

checkresiduals(USgas_arima_mod)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(2,1,1)[12]
## Q* = 24.949, df = 18, p-value = 0.1263
## 
## Model df: 6.   Total lags used: 24

Everything looks good here, we can see that the residuals resemble white noise and are evenly distributed with a mean of zero. There are no significant autocorrelations in the lags, and the Ljung-Box test also does not reject the null-hypothesis of no auto correlation with a non-significant p-value of 0.13.

Conclusion

And there we go, we’ve successfully fitted a few different models, both of which do a good job of capturing the data, but I think it’s fair to say that the ARIMA model has come out as the clear winner in this case. Let’s take another look at it’s predictions for the last 12 months, complete with a confidence interval:

plot_forecast(USgas_test_fc,
title = "US Natural Gas Consumption - Forecast",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")